% using Georg's approach/settings running onsets only for fb
%plots perturbation population response
%plots population resonse on running onset (fb)

ii=figure;
oo=figure;
perturb_act=[];
run_act=[];
l_win=10;
r_win=50;
for site=1:size(proj_meta,2)
no_time_points = size(proj_meta(site).rd,2);

act=[];
type_sessions={'f','p','d'};
for tp=1:no_time_points
    act=[];
    ps_onsets=[];
    ps_snip=[];
    cnt2=0;
    for lyr=1:size(proj_meta(site).rd,1)
        act=[act; proj_meta(site).rd(lyr,tp).act];
    end
    
    % indices for feedback session
    sess_start=proj_meta(site).rd(lyr,tp).nbr_frames;
    sess_start=cumsum(sess_start);
    sess_start=[1 sess_start+1];
    ac_sess=find(strcmp(proj_meta(site).rd(lyr,tp).session,type_sessions(1))==1);
    indices=[];
    for kk=1:length(ac_sess)
        indices=[indices sess_start(ac_sess(kk)):sess_start(ac_sess(kk)+1)-1];
    end
    avg_act=mean(act(:,indices));
    ps=proj_meta(site).rd(lyr,tp).ps_id(indices);
    run=proj_meta(site).rd(lyr,tp).velM_smoothed(indices);
    run=run-mean(run);
    run=run/max(run);
    run=run>0.05;
    
    ps_onsets=find(diff(round(ps))>0.05);
    ps_onsets(ps_onsets<=l_win)=[];
    ps_onsets(ps_onsets>=length(run)-r_win)=[];
    for ind=1:length(ps_onsets)
        cnt2=cnt2+1;
        ps_snip(:,cnt2)=avg_act(ps_onsets(ind)-l_win:ps_onsets(ind)+r_win);
    end
%     perturb_act(site,tp,1:61) = mean(ps_snip(90:90+60,:),2);
    perturb_act(site,tp,1:l_win+r_win+1) = mean(ps_snip,2);
    
% % % running onsets
% % % from playback session
% %     ac_sess=find(strcmp(proj_meta(site).rd(lyr,tp).session,type_sessions(2))==1);
% %     indices=[];
% %     for kk=1:length(ac_sess)
% %         indices=[indices sess_start(ac_sess(kk)):sess_start(ac_sess(kk)+1)-1];
% %     end
% %     avg_act=mean(act(:,indices));
% %     run=proj_meta(site).rd(lyr,tp).velM_smoothed(indices);
% %     run=run-mean(run);
% %     run=run/max(run);
% %     
% %     r_onsets=find(diff(run>0.05));
% %     r_onsets(logical([0 diff(r_onsets)<30]))=[];
% %     r_onsets=intersect(r_onsets,find(diff(run>0.05)==1));
% %     r_onsets(r_onsets<l_win)=[];
% %     r_onsets(r_onsets>length(run)-r_win)=[];
% %     
% %     for ind=1:length(r_onsets)
% %         cnt2=cnt2+1;
% %         run_snip(:,cnt2)=avg_act(r_onsets(ind)-l_win:r_onsets(ind)+r_win);
% %     end
% %     run_act(site,tp,1:l_win+r_win+1) = mean(run_snip,2);
    
perturb_act(site,tp,:)=perturb_act(site,tp,:)/mean(squeeze(perturb_act(site,tp,1:5)));
% run_act(site,tp,:)=run_act(site,tp,:)/mean(squeeze(run_act(site,tp,1:5)));

figure(ii);
colors=[0 0 1;0 0 1;1 0 0;1 0 0;0 0 0;0 0 0;0 0 0;0 0 0;0 0 0;0 0 0;];
% set(0,'DefaultAxesColorOrder',colors(1:no_time_points,:));
subplot(ceil(size(proj_meta,2)/2),2,site);
hold on
plot(squeeze(perturb_act(site,tp,:))','Color',colors(tp,:))
% plot(perturb_act');
title([regexprep(proj_meta(site).animal,'\_','\\_') ' - ' num2str(proj_meta(site).ExpGroup(1))]);

% % figure(oo);
% % % set(0,'DefaultAxesColorOrder',colors (1:no_time_points,:));
% % subplot(ceil(size(proj_meta,2)/2),2,site);
% % hold on
% % plot(squeeze(run_act(site,tp,:))','Color',colors(tp,:))
% % % plot(perturb_act');
% % title([regexprep(proj_meta(site).animal,'\_','\\_') ' - ' num2str(proj_meta(site).ExpGroup(1))]);
end
end

% average tp - mismatch
end_points = [5, 6, 9, 9, 9, 10];
figure;
for site = 1:size(proj_meta,2)
    subplot(ceil(size(proj_meta,2)/2),2,site);
    hold on
    plot(mean(squeeze(perturb_act(site,1:2,:)))','b')
    plot(mean(squeeze(perturb_act(site,3:4,:)))','r')
    if end_points(site) == 5
        plot(squeeze(perturb_act(site,5,:)),'k')
    else
        plot(mean(squeeze(perturb_act(site,5:end_points(site),:)),1),'k')
    end
    title([regexprep(proj_meta(site).animal,'\_','\\_') ' - ' num2str(proj_meta(site).ExpGroup(1))]);
end


% average tp - running onsets
end_points = [5, 6, 9, 9, 9, 10];
figure;
for site = 1:size(proj_meta,2) 
subplot(ceil(size(proj_meta,2)/2),2,site);
hold on
plot(mean(squeeze(run_act(site,1:2,:)))','b')
plot(mean(squeeze(run_act(site,3:4,:)))','r')
if end_points(site) == 5
    plot(squeeze(run_act(site,5,:)),'k')
else
plot(mean(squeeze(run_act(site,5:end_points(site),:)),1),'k')
end
title([regexprep(proj_meta(site).animal,'\_','\\_') ' - ' num2str(proj_meta(site).ExpGroup(1))]);
end



% % figure;
% % hold on,
% % % perturb_act=perturb_act([1:3,5:7],:);
% % % run_act=run_act([1:3,5:7],:);
% % % plot(perturb_act','r')
% % % plot(run_act','b')
% % % plot(mean(perturb_act,1),'r.')
% % % plot(mean(run_act,1),'b.')
% % % plot((mean(perturb_act,1)+std(perturb_act)/sqrt(6))','k+','linewidth',2)
% % % plot((mean(perturb_act,1)-std(perturb_act)/sqrt(6))','k+','linewidth',2)
% % % plot((mean(run_act,1)-std(run_act)/sqrt(6))','k+','linewidth',2)
% % % plot((mean(run_act,1)+std(run_act)/sqrt(6))','k+','linewidth',2)
% % 
% % aa=figure;hold on
% % ac=(mean(perturb_act)-1)*100;
% % ac_sem=((std(perturb_act)/sqrt(6)))*100;
% % ac2=(mean(run_act)-1)*100;
% % ac_sem2=(std(run_act)/sqrt(6))*100;
% % area([-10/10:1/10:50/10 50/10:-1/10:-10/10],[ac + ac_sem fliplr(ac-ac_sem)],'Facecolor',[0.9 0.9 0.9],'LineStyle','none');
% % area([-10/10:1/10:50/10 50/10:-1/10:-10/10],[ac2 + ac_sem2 fliplr(ac2-ac_sem2)],'Facecolor',[0.9 0.9 0.9],'LineStyle','none');
% % 
% % plot(-10/10:1/10:50/10,ac,'Color',[244/256 164/256 96/256],'Linewidth',3)
% % 
% % plot(-10/10:1/10:50/10,ac2,'Color',[0 0 0],'Linewidth',3)
% % xlim([-1.1 5.2])
% % ylim([-2 7])
% % xlabel('time [s]')
% % ylabel('\DeltaF/F [%]')
% % 
